home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / QUATERN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-02-07  |  11.8 KB  |  576 lines
  1. /*****************************************************************************
  2. *                quatern.c
  3. *
  4. *  This module implements Quaternion algebra julia fractals.
  5. *
  6. *  This file was written by Pascal Massimino.
  7. *  Revised and updated for POV-Ray 3.x by Tim Wegner
  8. *
  9. *  from Persistence of Vision(tm) Ray Tracer
  10. *  Copyright 1996 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file. If
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. #include "frame.h"
  28. #include "povray.h"
  29. #include "vector.h"
  30. #include "povproto.h"
  31. #include "fractal.h"
  32. #include "quatern.h"
  33. #include "spheres.h"
  34.  
  35.  
  36.  
  37. /*****************************************************************************
  38. * Local preprocessor defines
  39. ******************************************************************************/
  40.  
  41. #ifndef Fractal_Tolerance
  42. #define Fractal_Tolerance 1e-8
  43. #endif
  44.  
  45. #define Deriv_z2(n1,n2,n3,n4)               \
  46. {                                           \
  47.   tmp = (n1)*x - (n2)*y - (n3)*z - (n4)*w;  \
  48.   (n2) = (n1)*y + x*(n2);                   \
  49.   (n3) = (n1)*z + x*(n3);                   \
  50.   (n4) = (n1)*w + x*(n4);                   \
  51.   (n1) = tmp;                               \
  52. }
  53.  
  54. #define Deriv_z3(n1,n2,n3,n4)              \
  55. {                                          \
  56.   dtmp = 2.0*((n2)*y + (n3)*z + (n4)*w);   \
  57.   dtmp2 = 6.0*x*(n1) - dtmp;               \
  58.   (n1) = ( (n1)*x3 - x*dtmp )*3.0;         \
  59.   (n2) = (n2)*x4 + y*dtmp2;                \
  60.   (n3) = (n3)*x4 + z*dtmp2;                \
  61.   (n4) = (n4)*x4 + w*dtmp2;                \
  62. }
  63.  
  64.  
  65. /*****************************************************************************
  66. * Local typedefs
  67. ******************************************************************************/
  68.  
  69.  
  70.  
  71. /*****************************************************************************
  72. * Static functions
  73. ******************************************************************************/
  74.  
  75.  
  76.  
  77. /*****************************************************************************
  78. * Local variables
  79. ******************************************************************************/
  80.  
  81.  
  82.  
  83. /*****************************************************************************
  84. *
  85. * FUNCTION
  86. *
  87. * INPUT
  88. *
  89. * OUTPUT
  90. *
  91. * RETURNS
  92. *
  93. * AUTHOR
  94. *
  95. *   Pascal Massimino
  96. *
  97. * DESCRIPTION
  98. *
  99. *   -
  100. *
  101. * CHANGES
  102. *
  103. *   Dec 1994 : Creation.
  104. *
  105. ******************************************************************************/
  106.  
  107. int Iteration_z3(point, Julia)
  108. VECTOR point;
  109. FRACTAL * Julia;
  110. {
  111.   int i;
  112.   DBL x, y, z, w;
  113.   DBL d, x2, tmp;
  114.   DBL Exit_Value;
  115.  
  116.   Sx[0] = x = point[X];
  117.   Sy[0] = y = point[Y];
  118.   Sz[0] = z = point[Z];
  119.   Sw[0] = w = (Julia->SliceDist  
  120.                   - Julia->Slice[X]*x 
  121.                   - Julia->Slice[Y]*y 
  122.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  123.   
  124.   Exit_Value = Julia->Exit_Value;
  125.  
  126.   for (i = 1; i <= Julia->n; ++i)
  127.   {
  128.     d = y * y + z * z + w * w;
  129.  
  130.     x2 = x * x;
  131.  
  132.     if ((d + x2) > Exit_Value)
  133.     {
  134.       return (FALSE);
  135.     }
  136.  
  137.     tmp = 3.0 * x2 - d;
  138.  
  139.     Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
  140.     Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
  141.     Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
  142.     Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
  143.   }
  144.  
  145.   return (TRUE);
  146. }
  147.  
  148.  
  149.  
  150. /*****************************************************************************
  151. *
  152. * FUNCTION
  153. *
  154. * INPUT
  155. *
  156. * OUTPUT
  157. *
  158. * RETURNS
  159. *
  160. * AUTHOR
  161. *
  162. *   Pascal Massimino
  163. *
  164. * DESCRIPTION
  165. *
  166. *   -
  167. *
  168. * CHANGES
  169. *
  170. *   Dec 1994 : Creation.
  171. *
  172. ******************************************************************************/
  173.  
  174. int Iteration_Julia(point, Julia)
  175. VECTOR point;
  176. FRACTAL * Julia;
  177. {
  178.   int i;
  179.   DBL x, y, z, w;
  180.   DBL d, x2;
  181.   DBL Exit_Value;
  182.  
  183.   Sx[0] = x = point[X];
  184.   Sy[0] = y = point[Y];
  185.   Sz[0] = z = point[Z];
  186.   Sw[0] = w = (Julia->SliceDist  
  187.                   - Julia->Slice[X]*x 
  188.                   - Julia->Slice[Y]*y 
  189.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  190.  
  191.   Exit_Value = Julia->Exit_Value;
  192.  
  193.   for (i = 1; i <= Julia->n; ++i)
  194.   {
  195.     d = y * y + z * z + w * w;
  196.  
  197.     x2 = x * x;
  198.  
  199.     if ((d + x2) > Exit_Value)
  200.     {
  201.       return (FALSE);
  202.     }
  203.  
  204.     x *= 2.0;
  205.  
  206.     Sy[i] = y = x * y + Julia->Julia_Parm[Y];
  207.     Sz[i] = z = x * z + Julia->Julia_Parm[Z];
  208.     Sw[i] = w = x * w + Julia->Julia_Parm[T];
  209.     Sx[i] = x = x2 - d + Julia->Julia_Parm[X];;
  210.  
  211.   }
  212.  
  213.   return (TRUE);
  214. }
  215.  
  216.  
  217.  
  218. /*****************************************************************************
  219. *
  220. * FUNCTION
  221. *
  222. * INPUT
  223. *
  224. * OUTPUT
  225. *
  226. * RETURNS
  227. *
  228. * AUTHOR
  229. *
  230. *   Pascal Massimino
  231. *
  232. * DESCRIPTION
  233. *
  234. * D_Iteration puts in *Dist a lower bound for the distance from *point to the
  235. * set
  236. *
  237. * CHANGES
  238. *
  239. *   Dec 1994 : Creation.
  240. *
  241. ******************************************************************************/
  242.  
  243. /*----------- Distance estimator + iterations ------------*/
  244.  
  245. int D_Iteration_z3(point, Julia, Dist)
  246. VECTOR point;
  247. FRACTAL * Julia;
  248. DBL * Dist;
  249. {
  250.   int i, j;
  251.   DBL Norm, d;
  252.   DBL xx, yy, zz;
  253.   DBL x, y, z, w;
  254.   DBL tmp, x2;
  255.   DBL Exit_Value;
  256.   DBL Pow;
  257.  
  258.   x = Sx[0] = point[X];
  259.   y = Sy[0] = point[Y];
  260.   z = Sz[0] = point[Z];
  261.   w = Sw[0] = (Julia->SliceDist  
  262.                   - Julia->Slice[X]*x 
  263.                   - Julia->Slice[Y]*y 
  264.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  265.  
  266.   Exit_Value = Julia->Exit_Value;
  267.  
  268.   for (i = 1; i <= Julia->n; i++)
  269.   {
  270.     d = y * y + z * z + w * w;
  271.  
  272.     x2 = x * x;
  273.  
  274.     if ((Norm = d + x2) > Exit_Value)
  275.     {
  276.       /* Distance estimator */
  277.  
  278.       x = Sx[0];
  279.       y = Sy[0];
  280.       z = Sz[0];
  281.       w = Sw[0];
  282.  
  283.       Pow = 1.0 / 3.0;
  284.  
  285.       for (j = 1; j < i; ++j)
  286.       {
  287.         xx = x * Sx[j] - y * Sy[j] - z * Sz[j] - w * Sw[j];
  288.         yy = x * Sy[j] + y * Sx[j] - z * Sw[j] + w * Sz[j];
  289.         zz = x * Sz[j] + y * Sw[j] + z * Sx[j] - w * Sy[j];
  290.         w  = x * Sw[j] - y * Sz[j] + z * Sy[j] + w * Sx[j];
  291.  
  292.         x = xx;
  293.         y = yy;
  294.         z = zz;
  295.  
  296.         Pow /= 3.0;
  297.       }
  298.  
  299.       *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
  300.  
  301.       return (FALSE);
  302.     }
  303.  
  304.     tmp = 3.0 * x2 - d;
  305.  
  306.     Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
  307.     Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
  308.     Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
  309.     Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
  310.   }
  311.  
  312.   *Dist = Precision;
  313.  
  314.   return (TRUE);
  315. }
  316.  
  317.  
  318.  
  319. /*****************************************************************************
  320. *
  321. * FUNCTION
  322. *
  323. * INPUT
  324. *
  325. * OUTPUT
  326. *
  327. * RETURNS
  328. *
  329. * AUTHOR
  330. *
  331. *   Pascal Massimino
  332. *
  333. * DESCRIPTION
  334. *
  335. *   -
  336. *
  337. * CHANGES
  338. *
  339. *   Dec 1994 : Creation.
  340. *
  341. ******************************************************************************/
  342.  
  343. int D_Iteration_Julia(point, Julia, Dist)
  344. VECTOR point;
  345. FRACTAL * Julia;
  346. DBL * Dist;
  347. {
  348.   int i, j;
  349.   DBL Norm, d;
  350.   DBL Exit_Value;
  351.   DBL x, y, z, w;
  352.   DBL xx, yy, zz, x2;
  353.   DBL Pow;
  354.  
  355.   x = Sx[0] = point[X];
  356.   y = Sy[0] = point[Y];
  357.   z = Sz[0] = point[Z];
  358.   w = Sw[0] = (Julia->SliceDist  
  359.                   - Julia->Slice[X]*x 
  360.                   - Julia->Slice[Y]*y 
  361.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  362.  
  363.   Exit_Value = Julia->Exit_Value;
  364.  
  365.   for (i = 1; i <= Julia->n; i++)
  366.   {
  367.     d = y * y + z * z + w * w;
  368.  
  369.     x2 = x * x;
  370.  
  371.     if ((Norm = d + x2) > Exit_Value)
  372.     {
  373.       /* Distance estimator */
  374.  
  375.       x = Sx[0];
  376.       y = Sy[0];
  377.       z = Sz[0];
  378.       w = Sw[0];
  379.  
  380.       Pow = 1.0 / 2.0;
  381.  
  382.       for (j = 1; j < i; ++j)
  383.       {
  384.         xx = x * Sx[j] - y * Sy[j] - z * Sz[j] - w * Sw[j];
  385.         yy = x * Sy[j] + y * Sx[j] + w * Sz[j] - z * Sw[j];
  386.         zz = x * Sz[j] + z * Sx[j] + y * Sw[j] - w * Sy[j];
  387.         w  = x * Sw[j] + w * Sx[j] + z * Sy[j] - y * Sz[j];
  388.  
  389.         x = xx;
  390.         y = yy;
  391.         z = zz;
  392.  
  393.         Pow /= 2.0;
  394.       }
  395.  
  396.       *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
  397.  
  398.       return (FALSE);
  399.     }
  400.  
  401.     x *= 2.0;
  402.  
  403.     Sy[i] = y = x * y + Julia->Julia_Parm[Y];
  404.     Sz[i] = z = x * z + Julia->Julia_Parm[Z];
  405.     Sw[i] = w = x * w + Julia->Julia_Parm[T];
  406.     Sx[i] = x = x2 - d + Julia->Julia_Parm[X];
  407.  
  408.   }
  409.  
  410.   *Dist = Precision;
  411.  
  412.   return (TRUE);
  413. }
  414.  
  415. /*****************************************************************************
  416. *
  417. * FUNCTION
  418. *
  419. * INPUT
  420. *   
  421. * OUTPUT
  422. *   
  423. * RETURNS
  424. *
  425. * AUTHOR
  426. *
  427. *   Pascal Massimino
  428. *
  429. * DESCRIPTION
  430. *
  431. * Provided the iterations sequence has been built, perform the computation of
  432. * the Normal
  433. *
  434. * CHANGES
  435. *
  436. *   Dec 1994 : Creation.
  437. *
  438. ******************************************************************************/
  439.  
  440. void Normal_Calc_z3(Result, N_Max, fractal)
  441. VECTOR Result;
  442. int N_Max;
  443. FRACTAL *fractal;
  444. {
  445.   DBL
  446.   n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
  447.   n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
  448.   n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
  449.  
  450.   DBL x, y, z, w;
  451.   int i;
  452.   DBL tmp, dtmp, dtmp2, x2, x3, x4;
  453.  
  454.   x = Sx[0];
  455.   y = Sy[0];
  456.   z = Sz[0];
  457.   w = Sw[0];
  458.  
  459.   for (i = 1; i <= N_Max; i++)
  460.   {
  461.     tmp = y * y + z * z + w * w;
  462.  
  463.     x2 = x * x;
  464.     x3 = x2 - tmp;
  465.     x4 = 3.0 * x2 - tmp;
  466.  
  467.     Deriv_z3(n11, n12, n13, n14);
  468.     Deriv_z3(n21, n22, n23, n24);
  469.     Deriv_z3(n31, n32, n33, n34);
  470.  
  471.     x = Sx[i];
  472.     y = Sy[i];
  473.     z = Sz[i];
  474.     w = Sw[i];
  475.   }
  476.  
  477.   Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
  478.   Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
  479.   Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
  480. }
  481.  
  482.  
  483.  
  484. /*****************************************************************************
  485. *
  486. * FUNCTION
  487. *
  488. * INPUT
  489. *
  490. * OUTPUT
  491. *
  492. * RETURNS
  493. *
  494. * AUTHOR
  495. *
  496. *   Pascal Massimino
  497. *
  498. * DESCRIPTION
  499. *
  500. *   -
  501. *
  502. * CHANGES
  503. *
  504. *   Dec 1994 : Creation.
  505. *
  506. ******************************************************************************/
  507.  
  508. void Normal_Calc_Julia(Result, N_Max, fractal)
  509. VECTOR Result;
  510. int N_Max;
  511. FRACTAL *fractal;
  512. {
  513.   DBL
  514.   n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
  515.   n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
  516.   n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
  517.   DBL tmp;
  518.   DBL x, y, z, w;
  519.   int i;
  520.  
  521.   x = Sx[0];
  522.   y = Sy[0];
  523.   z = Sz[0];
  524.   w = Sw[0];
  525.  
  526.   for (i = 1; i <= N_Max; i++)
  527.   {
  528.     Deriv_z2(n11, n12, n13, n14);
  529.     Deriv_z2(n21, n22, n23, n24);
  530.     Deriv_z2(n31, n32, n33, n34);
  531.  
  532.     x = Sx[i];
  533.     y = Sy[i];
  534.     z = Sz[i];
  535.     w = Sw[i];
  536.   }
  537.  
  538.   Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
  539.   Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
  540.   Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
  541. }
  542.  
  543.  
  544.  
  545. /*****************************************************************************
  546. *
  547. * FUNCTION
  548. *
  549. * INPUT
  550. *
  551. * OUTPUT
  552. *   
  553. * RETURNS
  554. *   
  555. * AUTHOR
  556. *
  557. *   Pascal Massimino
  558. *   
  559. * DESCRIPTION
  560. *
  561. *   -
  562. *
  563. * CHANGES
  564. *
  565. *   Dec 1994 : Creation.
  566. *
  567. ******************************************************************************/
  568.  
  569. int F_Bound_Julia(Ray, Fractal, Depth_Min, Depth_Max)
  570. RAY *Ray;
  571. FRACTAL *Fractal;
  572. DBL *Depth_Min, *Depth_Max;
  573. {
  574.   return (Intersect_Sphere(Ray, Fractal->Center, Fractal->Radius_Squared, Depth_Min, Depth_Max));
  575. }
  576.